# A1 and A2
source('scripts/createDataset.R')

reg1 <- list()
reg2 <- list()
temp1 <- data.frame(stringsAsFactors = F)
balanceTestVar <- contVarBase
balanceTestVar_lab <- c("age", "married", '# children', 'w/o degree', 
                        "income (log)", 'work exp.', 'satisf. w/ HH-inc.')
for(i in 1:length(balanceTestVar)){
  pdf1 <- df %>% select(-y) %>% rename(y = balanceTestVar[i])
  contVarBase_temp <- contVarBase %>% subset(., !. %in% balanceTestVar[i])
  temp <- rdd_data(pdf1$y, pdf1$x, covar = pdf1[,contVarBase_temp], cutpoint = cutoff)
  reg1[[i]] <- temp %>% rdd_reg_lm(., slope = "separate", bw = bw, order = 1,
                                   covariates = contVarBase_temp %>% paste(., collapse = '+'))
  reg2[[i]] <- temp %>% rdd_reg_lm(., slope = "separate", bw = bw, order = 2,
                                   covariates = contVarBase_temp %>% paste(., collapse = '+'))
  temp1 <- temp %>%
    mutate(x = floor(x/6)*6) %>%
    group_by(x) %>% summarize(y = mean(y, na.rm = T), n = n()) %>% ungroup() %>%
    mutate(D = case_when(x < 0 ~ 'not eligible for Mütterrente', 
                         T ~ 'eligible for Mütterrente')) %>%
    mutate(var = balanceTestVar_lab[i]) %>%
    bind_rows(., temp1)
}

p <- ggplot() +
  geom_vline(xintercept = 0) +
  geom_point(aes(x = (-1)*(x+1), y = y, color = D), alpha = 0.9, stroke = 0, data = temp1) +
  geom_smooth(aes(x = (-1)*(x+1), y = y, color = D), size = 1.5, method = "lm", 
              formula = y ~ poly(x, 2), se = F, level = 0.95, alpha = 0.4, data = temp1) +
  scale_color_manual(values = c('eligible for Mütterrente' = 'black', 
                                'not eligible for Mütterrente' = 'gray50')) +
  scale_size_continuous(guide = "none") +
  scale_x_continuous(breaks = seq(-480, 480, 48)) +
  xlab('') + ylab('') +
  theme_light() + theme(legend.title = element_blank(), legend.position = 'none',
                        axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~var, scales = 'free')
ggsave("results/APPENDIX_A/A2.png", p, width = 7, height = 4)

fileName <- paste0("results/APPENDIX_A/A1.html")
names(reg2) <- balanceTestVar_lab
modelsummary(reg2, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

#############

# A3
source('scripts/createDataset.R')

x <- dat_temp_mccrary$x
rdd <- rddensity(X = x)
rdd[["test"]][["p_jk"]]
rdplotdensity(rdd, x, histFillCol = 'gray80', histFillShade = 1, CIshade = 0.4, lwd = 1) %>% .[["Estplot"]] +
  geom_vline(xintercept = 0) +
  scale_x_continuous(breaks = seq(-480, 480, 24)) +
  xlab('number of months prior/after January 1992') +
  ylab('density') +
  theme_bw() + 
  theme(legend.position = 'none')
ggsave("results/APPENDIX_A/A3.png", width = 8, height = 3)

#############

# A4
source('scripts/createDataset.R')

temp <- dat_temp %>%
  subset(., abs(x) <= bw) %>%
  mutate(x = floor(x/6)*6) %>%
  group_by(x) %>% summarize(y = mean(y, na.rm = T), n = n()) %>% ungroup() %>%
  mutate(D = case_when(x < 0 ~ 'not eligible for Mütterrente', T ~ 'eligible for Mütterrente'))
ggplot() +
  geom_vline(xintercept = 0) +
  geom_smooth(aes(x = (-1)*(x+1), y = y, color = D), size = 1, method = "lm",
              formula = y ~ poly(x, 2), se = F, level = 0.95, data = temp, alpha = 0.5) +
  geom_smooth(aes(x = (-1)*(x+1), y = y, color = D), size = 1, method = "lm", linetype = 'dashed',
              formula = y ~ poly(x, 1), se = F, level = 0.95, data = temp, alpha = 0.5) +
  geom_point(aes(x = (-1)*(x+1), y = y, color = D), size = 3, shape = 16, stroke = 0, data = temp) +
  scale_color_manual(values = c('eligible for Mütterrente' = 'black', 
                                'not eligible for Mütterrente' = 'gray50')) +
  scale_size_continuous(guide = "none") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 5L), breaks = seq(0, 1, 0.1), limits = c(0, 0.8)) +
  scale_x_continuous(breaks = seq(-480, 480, 12)) +
  xlab('number of months prior/after January 1992') +
  ylab('alignment with the \npledge-making party') +
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom')
ggsave("results/APPENDIX_A/A4.png", width = 10, height = 4)

#############

# A5
source('scripts/createDataset.R')

selectID <- df %>% 
  mutate(ID = 1:nrow(.)) %>%
  select(all_of(c(dep, run, contVarBase)), ID) %>%
  subset(., complete.cases(.)) %>%
  pull(ID)
temp <- df %>% 
  mutate(ID = 1:nrow(.)) %>%
  subset(., ID %in% selectID) %>%
  select(all_of(c(dep, 'x', contVarBase))) %>%
  as.data.frame() %>%
  mutate(SAMPLE = case_when(abs(x) <= 165 ~ 'in-sample', T ~ 'out-of-sample'))

t_test_results_df <- data.frame(stringsAsFactors = F)
for (var in contVarBase) {
  t_test_result <- t.test(as.formula(paste0(var, " ~ SAMPLE")), data = temp)
  means <- aggregate(as.formula(paste0(var, " ~ SAMPLE")), data = temp, FUN = mean)
  t_test_results_df <- rbind(t_test_results_df, data.frame(
    variable = var,
    insample = means[means$SAMPLE == "in-sample", 2],
    outsample = means[means$SAMPLE == "out-of-sample", 2],
    t = t_test_result$statistic,
    p = t_test_result$p.value
  ))
}
fileName <- paste0("results/APPENDIX_A/A5.html")
stargazer(t_test_results_df, type = "html", summary = FALSE, rownames = FALSE, out = fileName)

#############

# A6
source('scripts/createDataset.R')

df1 <- df
outcomes <- c("OTHgovtasksick", "OTHgovtaskfinsick", "OTHgovtaskunemployed", "OTHgovtaskjobs", 
              "OTHworrycrime", "OTHworryterrorism", "OTHworryeuro", 'OTHworryimmigration')
slope_temp <- c("separate")
order_temp <- c(1)

reg <- list()
for(p in 1:length(outcomes)){

  outcome_temp <- outcomes[p]
  df <- df1 %>% select(-y) %>% rename(y = outcome_temp)
  
  contVar <- c("age", "married", 'nchild', "nodegree", "income", "workexperience", "satishhinc")
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[length(reg)+1]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                                   covariates = contVar %>% paste(., collapse = '+'))
  reg[[length(reg)]]$outcome <- outcome_temp
}

fileName <- 'results/APPENDIX_A/A6.html'
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% rbind(c('outcome', unlist(lapply(reg, function(x) x$outcome)) %>% as.character()))
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, add_rows = info,
             notes = list(note2, note1), output = fileName)
